Regression (Phase-based) CLEAN

Intro

Here we provide exploration of LMER models fitted to the 2020 U.S. Congressional election Twitter Dataset (https://doi.org/10.6084/m9.figshare.25523734).

This document was typeset assembled using Quarto (https://quarto.org).

Initialize Workspace

Clear and Init Workspace

Imports

Code
library(ggplot2)
library(GGally)
library(plyr)
library(tidyverse)

Set Directories

Data Load and Inspection

Data Load

Code
ifn0 <- "i0001-politicians-aggregated.rds"
suppressWarnings(rm(list = ls(pattern = "^df")))
df0 <- readr::read_rds(file.path(ifd0, ifn0)) %>%
  dplyr::select(Name, Party, Outcome, Period, Phase, Time, Days, Agency) %>%
  dplyr::mutate(
    Period = forcats::fct_recode(
      Period,
      PC = "MI",
      PB = "BE",
      PA = "AR")) %>%
  identity()
df0 %>% dim() %>% cat0()
169997 8 
Code
for (name in names(df0)) {cat0(paste0("    ", name, ", "))}
    Name,  
    Party,  
    Outcome,  
    Period,  
    Phase,  
    Time,  
    Days,  
    Agency,  

Skim

Code
set.w(222)
library(skimr)
options(scipen=999)
skimr::skim_format(numeric=list(digits=6))
Warning: Formatting options are now in print() or as function arguments in skim_with().
Code
## digits=4, scientific = FALSE
df0 %>% skimr::skim() %>% dplyr::select(-c(n_missing)) %>% print(numeric=list(digits=2))
── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             169997    
Number of columns          8         
_______________________              
Column type frequency:               
  factor                   5         
  numeric                  3         
________________________             
Group variables            None      

── Variable type: factor ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate ordered n_unique top_counts                                
1 Name                      1 FALSE        870 Bar: 361, Ela: 361, Llo: 361, Pra: 361    
2 Party                     1 FALSE          2 Dem: 99081, Rep: 70916                    
3 Outcome                   1 FALSE          2 win: 115031, los: 54966                   
4 Period                    1 FALSE          3 PB: 98660, PA: 47974, PC: 23363           
5 Phase                     1 FALSE          4 BE: 98660, AR: 47974, AE: 12011, BR: 11352

── Variable type: numeric ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate     mean      sd      p0     p25     p50    p75   p100 hist 
1 Time                      1  -0.0680   0.572   -1     -0.55   -0.122  0.428   1    ▇▇▇▆▆
2 Days                      1 -12.2    103.    -180    -99     -22     77     180    ▇▇▇▆▆
3 Agency                    1   0.500    0.271   -1.69   0.343   0.500  0.658   2.37 ▁▁▇▂▁

Check factors order

Code
df2 <- df0
df3 <- df0
summary(df0$Period)
   PC    PB    PA 
23363 98660 47974 
Code
df2$Period <- factor(df2$Period, levels = c("PB", "PC", "PA"))
summary(df2$Period)
   PB    PC    PA 
98660 23363 47974 
Code
df3$Period <- relevel(df0$Period, ref = "PB")
summary(df3$Period)
   PB    PC    PA 
98660 23363 47974 

Phase Levels

Code
summary(df2$Phase) ##  <- factor(df2$Phase, levels = c("PB", "PC", "PA"))
   BE    AE    BR    AR 
98660 12011 11352 47974 
Code
attributes(df2$Phase)
$levels
[1] "BE" "AE" "BR" "AR"

$class
[1] "factor"

Check Phase and Period

Code
df0 %>%
  dplyr::group_by(Phase, Period) %>%
  dplyr::summarize(
    DaysMin = min(Days),
    DaysMax = max(Days),
    .groups = "drop")
# A tibble: 4 × 4
  Phase Period DaysMin DaysMax
  <fct> <fct>    <dbl>   <dbl>
1 BE    PB        -180       0
2 AE    PC           1      32
3 BR    PC          33      63
4 AR    PA          64     180

Prepare for Regression

File names and output directories

## ./data/j1002-model-by-period/i0001_models/fit0x0.extension
## ./data/j1008-model-by-phase/i0001_models/fit0x0.extension

Stuff for plotting

Marginalization

Checkups

Save data frame for reference

Code
file <- file.path(dir8, "data-df0.rds")
df0 %>% readr::write_rds(file)
cat0(file)
./data/j1008-model-by-phase/i0001_models/data-df0.rds 

Model 01: Null

Fit

Code
model <- "fit01aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (1 | Name) + 1,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit01aPh: [df0] Agency ~ (1 | Name) + 1
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (1 | Name) + 1
   Data: df0
Control: control

REML criterion at convergence: 26631.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.9514 -0.5639 -0.0023  0.5683  7.3986 

Random effects:
 Groups   Name        Variance Std.Dev.
 Name     (Intercept) 0.006209 0.0788  
 Residual             0.067519 0.2598  
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value            Pr(>|t|)    
(Intercept)   0.498618   0.002807 826.426886   177.7 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit01aPh: [df0] Agency ~ (1 | Name) + 1
# R2 for Mixed Models

  Conditional R2: 0.084
     Marginal R2: 0.000
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit01aPh: [df0] Agency ~ (1 | Name) + 1
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.084
  Unadjusted ICC: 0.084
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit01aPh: [df0] Agency ~ (1 | Name) + 1
# ICC by Group

Group |   ICC
-------------
Name  | 0.084

Random Effects

Code
extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 02: Time

Fit

Code
model <- "fit02aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit02aPh: [df0] Agency ~ (Time | Name) + Time
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time
   Data: df0
Control: control

REML criterion at convergence: 24941.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0353 -0.5639 -0.0047  0.5664  7.3362 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006283 0.07927      
          Time        0.004088 0.06394  0.18
 Residual             0.066384 0.25765      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value             Pr(>|t|)    
(Intercept)   0.498029   0.002848 816.443563 174.884 < 0.0000000000000002 ***
Time         -0.010230   0.002634 709.116999  -3.884             0.000112 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
Time 0.187 
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit02aPh: [df0] Agency ~ (Time | Name) + Time
# R2 for Mixed Models

  Conditional R2: 0.102
     Marginal R2: 0.000
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit02aPh: [df0] Agency ~ (Time | Name) + Time
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.102
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit02aPh: [df0] Agency ~ (Time | Name) + Time
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

Code
## pp(model);
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
tbl0
Characteristic Beta 95% CI1 p-value
Time -0.01 -0.02, -0.01 <0.001
Name.sd__(Intercept) 0.08

Name.cor__(Intercept).Time 0.18

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit02aPh: [df0] Agency ~ (Time | Name) + Time
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter   | Coefficient |       SE |         95% CI | t(169991) |      p
--------------------------------------------------------------------------
(Intercept) |        0.50 | 2.85e-03 | [ 0.49,  0.50] |    174.88 | < .001
Time        |       -0.01 | 2.63e-03 | [-0.02, -0.01] |     -3.88 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

Code
extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est") 
gg88 <- gg88 + ggplot2::ylim(-0.1,0.1) + line2 + stuff0
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

Code
terms <- c("Time")
extra <- "0001"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

pp(model);
fit02aPh: [df0] Agency ~ (Time | Name) + Time
Code
cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.51 | 0.50, 0.51
-0.50 |      0.50 | 0.50, 0.51
 0.00 |      0.50 | 0.49, 0.50
 0.50 |      0.50 | 0.49, 0.50
 1.00 |      0.49 | 0.48, 0.50
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-7.57e-03 | -0.01,  0.00 | 0.004
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.15,)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0
## gg88 <- gg88 + coord_cartesian(ylim = c(0.45, 0.55))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

Random Effects

Code
extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 03: Time x Phase

Fit

Code
model <- "fit03aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Phase,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Phase
   Data: df0
Control: control

REML criterion at convergence: 24187.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1124 -0.5632 -0.0069  0.5646  7.3739 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006236 0.07897      
          Time        0.004124 0.06422  0.18
 Residual             0.066073 0.25705      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                  Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)       0.528553      0.003172   1283.400354 166.622 < 0.0000000000000002 ***
Time              0.044610      0.003657   2775.135404  12.198 < 0.0000000000000002 ***
PhaseAE           0.015184      0.004861 168870.559183   3.123              0.00179 ** 
PhaseBR          -0.132158      0.012989 168617.283756 -10.175 < 0.0000000000000002 ***
PhaseAR          -0.000602      0.004747 169706.699223  -0.127              0.89909    
Time:PhaseAE     -0.638482      0.044797 168682.837032 -14.253 < 0.0000000000000002 ***
Time:PhaseBR      0.224563      0.047871 168575.449036   4.691           0.00000272 ***
Time:PhaseAR     -0.110523      0.007010 167010.782296 -15.765 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PhasAE PhasBR PhasAR Tm:PAE Tm:PBR
Time         0.431                                          
PhaseAE     -0.165 -0.219                                   
PhaseBR     -0.061 -0.083  0.043                            
PhaseAR     -0.171 -0.223  0.117  0.044                     
Time:PhasAE -0.026 -0.049 -0.786  0.009  0.022              
Time:PhasBR -0.026 -0.045  0.017 -0.966  0.019  0.004       
Time:PhasAR -0.172 -0.314  0.117  0.046 -0.695  0.030  0.025
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
# R2 for Mixed Models

  Conditional R2: 0.106
     Marginal R2: 0.005
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.101
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Code
tbl0
Characteristic Beta 95% CI1 p-value
Time 0.04 0.04, 0.05 <0.001
Phase


    AE - BE 0.06 0.04, 0.08 <0.001
    BR - BE -0.15 -0.19, -0.11 <0.001
    BR - AE -0.21 -0.25, -0.16 <0.001
    AR - BE 0.01 -0.01, 0.02 0.5
    AR - AE -0.05 -0.07, -0.03 <0.001
    AR - BR 0.15 0.11, 0.20 <0.001
Time * Phase


    Time * AE -0.64 -0.73, -0.55 <0.001
    Time * BR 0.22 0.13, 0.32 <0.001
    Time * AR -0.11 -0.12, -0.10 <0.001
Name.sd__(Intercept) 0.08

Name.cor__(Intercept).Time 0.18

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter         | Coefficient |       SE |         95% CI | t(169985) |      p
--------------------------------------------------------------------------------
(Intercept)       |        0.53 | 3.17e-03 | [ 0.52,  0.53] |    166.62 | < .001
Time              |        0.04 | 3.66e-03 | [ 0.04,  0.05] |     12.20 | < .001
Phase [AE]        |        0.02 | 4.86e-03 | [ 0.01,  0.02] |      3.12 | 0.002 
Phase [BR]        |       -0.13 |     0.01 | [-0.16, -0.11] |    -10.17 | < .001
Phase [AR]        |   -6.02e-04 | 4.75e-03 | [-0.01,  0.01] |     -0.13 | 0.899 
Time × Phase [AE] |       -0.64 |     0.04 | [-0.73, -0.55] |    -14.25 | < .001
Time × Phase [BR] |        0.22 |     0.05 | [ 0.13,  0.32] |      4.69 | < .001
Time × Phase [AR] |       -0.11 | 7.01e-03 | [-0.12, -0.10] |    -15.77 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

Code
pp(model)
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est")
gg88 <- gg88 + ggplot2::ylim(-0.3,0.3) + line2 + stuff0
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
Code
gg88
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

FE: Time

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "0002"
terms=c("Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "0003"
terms=c("Time", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

Code
pp(model)
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.52, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.52 | 0.52, 0.53
 0.50 |      0.52 | 0.51, 0.52
 1.00 |      0.51 | 0.50, 0.52
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.01 | -0.02,  0.00 | 0.025
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,) 
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Phase

NULL: Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "1002"
terms <- c("Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Phase | Predicted |     95% CI
------------------------------
BE    |      0.53 | 0.52, 0.53
AE    |      0.59 | 0.57, 0.60
BR    |      0.38 | 0.35, 0.41
AR    |      0.54 | 0.52, 0.55
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Phase | Predicted |     95% CI |      p
---------------------------------------
BE    |      0.53 | 0.52, 0.53 | < .001
AE    |      0.59 | 0.57, 0.60 | < .001
BR    |      0.38 | 0.35, 0.41 | < .001
AR    |      0.54 | 0.52, 0.55 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
## gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Phase |  Contrast |       95% CI |      p
-----------------------------------------
BE-AE |     -0.06 | -0.07, -0.04 | < .001
BE-BR |      0.15 |  0.12,  0.18 | < .001
BE-AR | -6.92e-03 | -0.02,  0.00 | 0.174 
AE-BR |      0.21 |  0.17,  0.24 | < .001
AE-AR |      0.05 |  0.03,  0.07 | < .001
BR-AR |     -0.15 | -0.19, -0.12 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase

NULL: Time x Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
extra <- "1003"
terms <- c("Time", "Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")


cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Phase: BE

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.47, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.56 | 0.55, 0.56
 1.00 |      0.58 | 0.57, 0.59

Phase: AE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      1.14 |  1.04, 1.23
-0.50 |      0.84 |  0.79, 0.89
 0.00 |      0.55 |  0.54, 0.56
 0.50 |      0.25 |  0.21, 0.29
 1.00 |     -0.04 | -0.12, 0.04

Phase: BR

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.13 | 0.01, 0.24
-0.50 |      0.26 | 0.19, 0.33
 0.00 |      0.40 | 0.37, 0.42
 0.50 |      0.54 | 0.51, 0.56
 1.00 |      0.67 | 0.60, 0.74

Phase: AR

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.59 | 0.57, 0.61
-0.50 |      0.56 | 0.55, 0.58
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.50 | 0.49, 0.51
 1.00 |      0.47 | 0.46, 0.48
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Phase | Slope |       95% CI |      p
-------------------------------------
BE    |  0.04 |  0.04,  0.05 | < .001
AE    | -0.59 | -0.68, -0.51 | < .001
BR    |  0.27 |  0.18,  0.36 | < .001
AR    | -0.07 | -0.08, -0.05 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
gg88 <- gg88 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Phase

Code
pp(model);
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |     0.64 |  0.55,  0.73 | < .001
BE-BR |    -0.22 | -0.32, -0.13 | < .001
BE-AR |     0.11 |  0.10,  0.12 | < .001
AE-BR |    -0.86 | -0.99, -0.73 | < .001
AE-AR |    -0.53 | -0.62, -0.44 | < .001
BR-AR |     0.34 |  0.24,  0.43 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

Model 04: Time x Phase x Outcome

Fit

Code
model <- "fit04aPh"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Phase * Outcome,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

Code
readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Phase * Outcome
   Data: df0
Control: control

REML criterion at convergence: 23685.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1858 -0.5634 -0.0081  0.5638  7.3838 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Name     (Intercept) 0.004804 0.06931       
          Time        0.003343 0.05782  -0.06
 Residual             0.065972 0.25685       
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                                Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)                     0.515092      0.004452   1470.335766 115.699 < 0.0000000000000002 ***
Time                            0.057026      0.005598   3432.607484  10.188 < 0.0000000000000002 ***
PhaseAE                        -0.017658      0.008833 169175.854063  -1.999              0.04559 *  
PhaseBR                        -0.145312      0.026930 168955.292605  -5.396         0.0000000683 ***
PhaseAR                        -0.104718      0.009831 166212.881574 -10.652 < 0.0000000000000002 ***
Outcomewinner                   0.023922      0.005856   1424.761477   4.085         0.0000465435 ***
Time:PhaseAE                   -0.882126      0.086357 169124.880897 -10.215 < 0.0000000000000002 ***
Time:PhaseBR                    0.017444      0.099141 168923.961312   0.176              0.86033    
Time:PhaseAR                   -0.077030      0.014696 150506.476445  -5.241         0.0000001596 ***
Time:Outcomewinner             -0.020359      0.007198   3282.626151  -2.828              0.00471 ** 
PhaseAE:Outcomewinner           0.051078      0.010583 169008.316334   4.826         0.0000013927 ***
PhaseBR:Outcomewinner           0.022108      0.030740 168862.263648   0.719              0.47202    
PhaseAR:Outcomewinner           0.140686      0.011241 167976.986518  12.516 < 0.0000000000000002 ***
Time:PhaseAE:Outcomewinner      0.306640      0.101049 168997.306098   3.035              0.00241 ** 
Time:PhaseBR:Outcomewinner      0.281534      0.113190 168837.528411   2.487              0.01287 *  
Time:PhaseAR:Outcomewinner     -0.035756      0.016747 157881.162351  -2.135              0.03276 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of fixed effects could have been required in summary()

Correlation of Fixed Effects:
            (Intr) Time   PhasAE PhasBR PhasAR Otcmwn Tm:PAE Tm:PBR Tm:PAR Tm:Otc PhAE:O PhBR:O PhAR:O T:PAE: T:PBR:
Time         0.364                                                                                                  
PhaseAE     -0.156 -0.194                                                                                           
PhaseBR     -0.051 -0.064  0.029                                                                                    
PhaseAR     -0.141 -0.177  0.079  0.026                                                                             
Outcomewnnr -0.760 -0.277  0.119  0.039  0.108                                                                      
Time:PhasAE -0.024 -0.043 -0.783  0.008  0.017  0.018                                                               
Time:PhasBR -0.021 -0.037  0.011 -0.973  0.013  0.016  0.002                                                        
Time:PhasAR -0.145 -0.248  0.075  0.027 -0.794  0.110  0.023  0.016                                                 
Tm:Otcmwnnr -0.283 -0.778  0.151  0.050  0.138  0.355  0.033  0.029  0.193                                          
PhsAE:Otcmw  0.130  0.162 -0.835 -0.024 -0.066 -0.169  0.653 -0.009 -0.063 -0.211                                   
PhsBR:Otcmw  0.045  0.056 -0.025 -0.876 -0.023 -0.058 -0.007  0.852 -0.024 -0.073  0.034                            
PhsAR:Otcmw  0.124  0.155 -0.069 -0.023 -0.875 -0.161 -0.015 -0.012  0.694 -0.200  0.094  0.032                     
Tm:PhsAE:Ot  0.020  0.037  0.669 -0.007 -0.014 -0.026 -0.855 -0.002 -0.020 -0.047 -0.783  0.008  0.018              
Tm:PhsBR:Ot  0.019  0.032 -0.010  0.852 -0.012 -0.024 -0.002 -0.876 -0.014 -0.041  0.014 -0.970  0.015  0.003       
Tm:PhsAR:Ot  0.127  0.218 -0.066 -0.024  0.696 -0.162 -0.021 -0.014 -0.878 -0.280  0.092  0.034 -0.759  0.025  0.019
Code
cat0(sep2); pp(model); performance::r2(get(model))
--------------------------------------------------------------------- 
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
# R2 for Mixed Models

  Conditional R2: 0.101
     Marginal R2: 0.020
Code
cat0(sep2); pp(model); performance::icc(get(model))
--------------------------------------------------------------------- 
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.083
  Unadjusted ICC: 0.081
Code
cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)
--------------------------------------------------------------------- 
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Model contains random slopes. Cannot compute accurate ICCs by group factors.
# ICC by Group

Group |   ICC
-------------
Name  | 0.067

Model Summary

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.
NOTE: Results may be misleading due to involvement in interactions
Code
tbl0
Characteristic Beta 95% CI1 p-value
Time 0.06 0.05, 0.07 <0.001
Phase


    AE - BE 0.06 0.04, 0.08 <0.001
    BR - BE -0.15 -0.19, -0.10 <0.001
    BR - AE -0.20 -0.26, -0.15 <0.001
    AR - BE -0.03 -0.04, -0.01 <0.001
    AR - AE -0.09 -0.11, -0.06 <0.001
    AR - BR 0.12 0.07, 0.17 <0.001
Outcome


    winner - loser 0.07 0.05, 0.09 <0.001
Time * Phase


    Time * AE -0.88 -1.1, -0.71 <0.001
    Time * BR 0.02 -0.18, 0.21 0.9
    Time * AR -0.08 -0.11, -0.05 <0.001
Time * Outcome


    Time * winner -0.02 -0.03, -0.01 0.005
Phase * Outcome


    AE * winner 0.05 0.03, 0.07 <0.001
    BR * winner 0.02 -0.04, 0.08 0.5
    AR * winner 0.14 0.12, 0.16 <0.001
Time * Phase * Outcome


    Time * AE * winner 0.31 0.11, 0.50 0.002
    Time * BR * winner 0.28 0.06, 0.50 0.013
    Time * AR * winner -0.04 -0.07, 0.00 0.033
Name.sd__(Intercept) 0.07

Name.cor__(Intercept).Time -0.06

Name.sd__Time 0.06

Residual.sd__Observation 0.26

1 CI = Confidence Interval

Model Parameters

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
parameters::model_parameters(get(model))
# Fixed Effects

Parameter                              | Coefficient |       SE |         95% CI | t(169977) |      p
-----------------------------------------------------------------------------------------------------
(Intercept)                            |        0.52 | 4.45e-03 | [ 0.51,  0.52] |    115.70 | < .001
Time                                   |        0.06 | 5.60e-03 | [ 0.05,  0.07] |     10.19 | < .001
Phase [AE]                             |       -0.02 | 8.83e-03 | [-0.03,  0.00] |     -2.00 | 0.046 
Phase [BR]                             |       -0.15 |     0.03 | [-0.20, -0.09] |     -5.40 | < .001
Phase [AR]                             |       -0.10 | 9.83e-03 | [-0.12, -0.09] |    -10.65 | < .001
Outcome [winner]                       |        0.02 | 5.86e-03 | [ 0.01,  0.04] |      4.08 | < .001
Time × Phase [AE]                      |       -0.88 |     0.09 | [-1.05, -0.71] |    -10.21 | < .001
Time × Phase [BR]                      |        0.02 |     0.10 | [-0.18,  0.21] |      0.18 | 0.860 
Time × Phase [AR]                      |       -0.08 |     0.01 | [-0.11, -0.05] |     -5.24 | < .001
Time × Outcome [winner]                |       -0.02 | 7.20e-03 | [-0.03, -0.01] |     -2.83 | 0.005 
Phase [AE] × Outcome [winner]          |        0.05 |     0.01 | [ 0.03,  0.07] |      4.83 | < .001
Phase [BR] × Outcome [winner]          |        0.02 |     0.03 | [-0.04,  0.08] |      0.72 | 0.472 
Phase [AR] × Outcome [winner]          |        0.14 |     0.01 | [ 0.12,  0.16] |     12.52 | < .001
(Time × Phase [AE]) × Outcome [winner] |        0.31 |     0.10 | [ 0.11,  0.50] |      3.03 | 0.002 
(Time × Phase [BR]) × Outcome [winner] |        0.28 |     0.11 | [ 0.06,  0.50] |      2.49 | 0.013 
(Time × Phase [AR]) × Outcome [winner] |       -0.04 |     0.02 | [-0.07,  0.00] |     -2.14 | 0.033 

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.07
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |       -0.06
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

FE: Time

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0002"
terms=c("Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0


file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0003"
terms=c("Time", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0004"
terms=c("Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0005"
terms=c("Time", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + lineR + line1 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Phase x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0006"
terms=c("Phase", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome x Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "0007"
terms=c("Time", "Outcome", "Phase")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 


gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)
## gg88 <- gg88 + ggplot2::coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time:

NULL: Time

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.52, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.52 | 0.52, 0.53
 0.50 |      0.52 | 0.51, 0.52
 1.00 |      0.51 | 0.50, 0.52
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.02 | -0.03,  0.00 | 0.006
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(-0.25, 1.25))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Phase

NULL: Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1002"
terms <- c("Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Phase | Predicted |     95% CI
------------------------------
BE    |      0.52 | 0.52, 0.53
AE    |      0.60 | 0.58, 0.62
BR    |      0.39 | 0.36, 0.43
AR    |      0.52 | 0.51, 0.53
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Phase | Predicted |     95% CI |      p
---------------------------------------
BE    |      0.52 | 0.52, 0.53 | < .001
AE    |      0.60 | 0.58, 0.62 | < .001
BR    |      0.39 | 0.36, 0.43 | < .001
AR    |      0.52 | 0.51, 0.53 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(
  limit_range=FALSE,  
  show_data = FALSE, 
  dot_alpha = 0.05)

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |    -0.08 | -0.09, -0.06 | < .001
BE-BR |     0.13 |  0.10,  0.17 | < .001
BE-AR | 4.18e-03 | -0.01,  0.01 | 0.442 
AE-BR |     0.21 |  0.17,  0.25 | < .001
AE-AR |     0.08 |  0.06,  0.10 | < .001
BR-AR |    -0.13 | -0.17, -0.09 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase

NULL: Time x Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1003"
terms <- c("Phase", "Time")
terms <- c("Time", "Phase")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.48 | 0.48,  0.49
-0.50 |      0.51 | 0.50,  0.51
 0.00 |      0.53 | 0.52,  0.53
 0.50 |      0.55 | 0.54,  0.56
 1.00 |      0.57 | 0.56,  0.58

Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.17 |  1.08,  1.27
-0.50 |      0.86 |  0.81,  0.91
 0.00 |      0.55 |  0.54,  0.56
 0.50 |      0.23 |  0.19,  0.27
 1.00 |     -0.08 | -0.16,  0.00

Phase: BR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.14 | 0.02,  0.27
-0.50 |      0.27 | 0.20,  0.35
 0.00 |      0.40 | 0.37,  0.42
 0.50 |      0.53 | 0.50,  0.55
 1.00 |      0.65 | 0.58,  0.72

Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.57 | 0.55,  0.60
-0.50 |      0.55 | 0.53,  0.56
 0.00 |      0.52 | 0.51,  0.53
 0.50 |      0.49 | 0.49,  0.50
 1.00 |      0.46 | 0.46,  0.47
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Phase | Slope |       95% CI |      p
-------------------------------------
BE    |  0.06 |  0.05,  0.07 | < .001
AE    | -0.83 | -0.99, -0.66 | < .001
BR    |  0.07 | -0.12,  0.27 | 0.452 
AR    | -0.02 | -0.05,  0.01 | 0.219 
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Phase

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Time x Phase

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Phase | Contrast |       95% CI |      p
----------------------------------------
BE-AE |     0.88 |  0.71,  1.05 | < .001
BE-BR |    -0.02 | -0.21,  0.18 | 0.860 
BE-AR |     0.08 |  0.05,  0.11 | < .001
AE-BR |    -0.90 | -1.16, -0.64 | < .001
AE-AR |    -0.81 | -0.98, -0.63 | < .001
BR-AR |     0.09 | -0.10,  0.29 | 0.414 
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Outcome

NULL: Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1004"
terms <- c("Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome | Predicted |     95% CI
--------------------------------
loser   |      0.45 | 0.44, 0.46
winner  |      0.52 | 0.51, 0.52
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Outcome | Predicted |     95% CI |      p
-----------------------------------------
loser   |      0.45 | 0.44, 0.46 | < .001
winner  |      0.52 | 0.51, 0.52 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.07 | -0.08, -0.06 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Outcome:

NULL: Time x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1005"
terms <- c("Time", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.50 | 0.47, 0.52
-0.50 |      0.48 | 0.47, 0.50
 0.00 |      0.47 | 0.46, 0.48
 0.50 |      0.46 | 0.45, 0.47
 1.00 |      0.45 | 0.43, 0.47

Outcome: winner

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.55 | 0.54, 0.57
-0.50 |      0.55 | 0.54, 0.56
 0.00 |      0.54 | 0.53, 0.55
 0.50 |      0.53 | 0.53, 0.54
 1.00 |      0.53 | 0.51, 0.54
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Outcome | Slope |     95% CI |      p
-------------------------------------
loser   |  0.06 | 0.05, 0.07 | < .001
winner  |  0.04 | 0.03, 0.05 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT Time x Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Outcome      | Contrast |     95% CI |     p
--------------------------------------------
loser-winner |     0.02 | 0.01, 0.03 | 0.005
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Phase x Outcome

NULL: Phase x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1006"
terms <- c("Phase", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser

Phase | Predicted |     95% CI
------------------------------
BE    |      0.51 | 0.50, 0.52
AE    |      0.55 | 0.52, 0.58
BR    |      0.36 | 0.30, 0.43
AR    |      0.41 | 0.39, 0.43

Outcome: winner

Phase | Predicted |     95% CI
------------------------------
BE    |      0.53 | 0.53, 0.54
AE    |      0.61 | 0.59, 0.62
BR    |      0.39 | 0.35, 0.43
AR    |      0.58 | 0.57, 0.59
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
Phase | Outcome | Predicted |     95% CI |      p
-------------------------------------------------
BE    |   loser |      0.51 | 0.50, 0.52 | < .001
AE    |   loser |      0.55 | 0.52, 0.58 | < .001
BR    |   loser |      0.36 | 0.30, 0.43 | < .001
AR    |   loser |      0.41 | 0.39, 0.43 | < .001
BE    |  winner |      0.53 | 0.53, 0.54 | < .001
AE    |  winner |      0.61 | 0.59, 0.62 | < .001
BR    |  winner |      0.39 | 0.35, 0.43 | < .001
AR    |  winner |      0.58 | 0.57, 0.59 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Phase x Outcome

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

PLOT: Phase x Outcome [low-level]

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1006"
terms <- c("Phase", "Outcome")

type <- "eff"
type <- "pred"

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- sjPlot::plot_model(
  get(model), 
  type = type, 
  terms = terms, 
  )

gg44 <- gg44 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg44, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg44, width=8, height=6)
gg44

DIFF: Phase x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# Pairwise comparisons

Phase |       Outcome | Contrast |       95% CI |      p
--------------------------------------------------------
BE-AE |   loser-loser |    -0.04 | -0.07, -0.02 | 0.003 
BE-BR |   loser-loser |     0.15 |  0.08,  0.21 | < .001
BE-AR |   loser-loser |     0.10 |  0.08,  0.12 | < .001
BE-BE |  loser-winner |    -0.03 | -0.04, -0.01 | < .001
BE-AE |  loser-winner |    -0.10 | -0.12, -0.08 | < .001
BE-BR |  loser-winner |     0.12 |  0.08,  0.16 | < .001
BE-AR |  loser-winner |    -0.07 | -0.08, -0.05 | < .001
AE-BR |   loser-loser |     0.19 |  0.12,  0.26 | < .001
AE-AR |   loser-loser |     0.14 |  0.11,  0.18 | < .001
AE-BE |  loser-winner |     0.02 | -0.01,  0.05 | 0.265 
AE-AE |  loser-winner |    -0.06 | -0.09, -0.02 | 0.001 
AE-BR |  loser-winner |     0.16 |  0.11,  0.21 | < .001
AE-AR |  loser-winner |    -0.03 | -0.06,  0.00 | 0.103 
BR-AR |   loser-loser |    -0.05 | -0.12,  0.02 | 0.201 
BR-BE |  loser-winner |    -0.17 | -0.24, -0.11 | < .001
BR-AE |  loser-winner |    -0.24 | -0.31, -0.18 | < .001
BR-BR |  loser-winner |    -0.03 | -0.10,  0.05 | 0.463 
BR-AR |  loser-winner |    -0.22 | -0.28, -0.15 | < .001
AR-BE |  loser-winner |    -0.12 | -0.15, -0.10 | < .001
AR-AE |  loser-winner |    -0.20 | -0.23, -0.17 | < .001
AR-BR |  loser-winner |     0.02 | -0.02,  0.06 | 0.400 
AR-AR |  loser-winner |    -0.17 | -0.19, -0.14 | < .001
BE-AE | winner-winner |    -0.07 | -0.09, -0.06 | < .001
BE-BR | winner-winner |     0.14 |  0.11,  0.18 | < .001
BE-AR | winner-winner |    -0.04 | -0.06, -0.03 | < .001
AE-BR | winner-winner |     0.22 |  0.18,  0.26 | < .001
AE-AR | winner-winner |     0.03 |  0.01,  0.05 | 0.006 
BR-AR | winner-winner |    -0.19 | -0.22, -0.15 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Phase x Outcome

NULL: Time x Phase x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
extra <- "1007"
terms <- c("Time", "Phase", "Outcome")
terms <- c("Time", "Outcome", "Phase") ## CAUTION
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)
===================================================================== 
Code
pred0 %>% print(n=Inf)
# Average predicted values of Agency

Outcome: loser
Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.45 | 0.44,  0.46
-0.50 |      0.48 | 0.47,  0.49
 0.00 |      0.51 | 0.50,  0.52
 0.50 |      0.54 | 0.53,  0.55
 1.00 |      0.57 | 0.56,  0.59

Outcome: loser
Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.32 |  1.13,  1.50
-0.50 |      0.91 |  0.81,  1.01
 0.00 |      0.49 |  0.48,  0.51
 0.50 |      0.08 |  0.01,  0.16
 1.00 |     -0.33 | -0.48, -0.17

Outcome: loser
Phase: BR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.29 | 0.04,  0.54
-0.50 |      0.33 | 0.18,  0.48
 0.00 |      0.37 | 0.31,  0.42
 0.50 |      0.41 | 0.36,  0.45
 1.00 |      0.44 | 0.30,  0.59

Outcome: loser
Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.43 | 0.38,  0.47
-0.50 |      0.42 | 0.38,  0.45
 0.00 |      0.41 | 0.39,  0.43
 0.50 |      0.40 | 0.39,  0.41
 1.00 |      0.39 | 0.37,  0.41

Outcome: winner
Phase: BE

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.50 | 0.49,  0.51
-0.50 |      0.52 | 0.51,  0.52
 0.00 |      0.54 | 0.53,  0.54
 0.50 |      0.56 | 0.55,  0.57
 1.00 |      0.58 | 0.56,  0.59

Outcome: winner
Phase: AE

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      1.11 |  0.99,  1.22
-0.50 |      0.84 |  0.78,  0.90
 0.00 |      0.57 |  0.56,  0.58
 0.50 |      0.30 |  0.26,  0.34
 1.00 |      0.03 | -0.06,  0.13

Outcome: winner
Phase: BR

 Time | Predicted |       95% CI
--------------------------------
-1.00 |      0.08 | -0.06,  0.21
-0.50 |      0.24 |  0.16,  0.33
 0.00 |      0.41 |  0.38,  0.44
 0.50 |      0.58 |  0.56,  0.61
 1.00 |      0.75 |  0.67,  0.83

Outcome: winner
Phase: AR

 Time | Predicted |      95% CI
-------------------------------
-1.00 |      0.65 | 0.62,  0.67
-0.50 |      0.61 | 0.59,  0.63
 0.00 |      0.57 | 0.56,  0.58
 0.50 |      0.54 | 0.53,  0.54
 1.00 |      0.50 | 0.49,  0.51
Code
cat0(sep0)
===================================================================== 
Code
test0 %>% print(n=Inf)
# (Average) Linear trend for Time

Outcome | Phase | Slope |       95% CI |      p
-----------------------------------------------
loser   |    BE |  0.06 |  0.05,  0.07 | < .001
loser   |    AE | -0.83 | -0.99, -0.66 | < .001
loser   |    BR |  0.07 | -0.12,  0.27 | 0.452 
loser   |    AR | -0.02 | -0.05,  0.01 | 0.187 
winner  |    BE |  0.04 |  0.03,  0.05 | < .001
winner  |    AE | -0.54 | -0.64, -0.44 | < .001
winner  |    BR |  0.34 |  0.23,  0.44 | < .001
winner  |    AR | -0.08 | -0.09, -0.06 | < .001
Code
pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Phase x Outcome

Code
extra <- "1007"
terms <- c("Time", "Phase", "Outcome")
terms <- c("Time", "Outcome", "Phase") ## CAUTION
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 

gg88 <- gg88 + rect0

# gg88 <- gg88 + ggplot2::ylim(0.10,0.90)
gg88 <- gg88 + scale_y_continuous(breaks = seq(0,1,by=0.1))
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Code
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=16)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=16)
gg88

Code
knitr::opts_chunk$set(fig.width=unit(12,"cm"), fig.height=unit(18,"cm"))

gg88

PLOT: Time x Phase x Outcome (low-level)

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
          get(model),
          type="emm",
          terms=terms)

gg88 <- gg88 + line0 + lineR + line1 + time0 + stuff0 + rect0
## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=4, height=16)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=4, height=16)
gg88

DIFF: Time x Phase x Outcome

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()
# (Average) Linear trend for Time

Outcome       | Phase | Contrast |       95% CI |      p
--------------------------------------------------------
loser-loser   | BE-AE |     0.88 |  0.71,  1.05 | < .001
loser-loser   | BE-BR |    -0.02 | -0.21,  0.18 | 0.860 
loser-loser   | BE-AR |     0.08 |  0.05,  0.11 | < .001
loser-winner  | BE-BE |     0.02 |  0.01,  0.03 | 0.006 
loser-winner  | BE-AE |     0.60 |  0.49,  0.70 | < .001
loser-winner  | BE-BR |    -0.28 | -0.39, -0.17 | < .001
loser-winner  | BE-AR |     0.13 |  0.11,  0.15 | < .001
loser-loser   | AE-BR |    -0.90 | -1.16, -0.64 | < .001
loser-loser   | AE-AR |    -0.81 | -0.98, -0.63 | < .001
loser-winner  | AE-BE |    -0.86 | -1.03, -0.69 | < .001
loser-winner  | AE-AE |    -0.29 | -0.48, -0.09 | 0.006 
loser-winner  | AE-BR |    -1.16 | -1.36, -0.96 | < .001
loser-winner  | AE-AR |    -0.75 | -0.92, -0.58 | < .001
loser-loser   | BR-AR |     0.09 | -0.10,  0.29 | 0.371 
loser-winner  | BR-BE |     0.04 | -0.16,  0.23 | 0.729 
loser-winner  | BR-AE |     0.61 |  0.39,  0.83 | < .001
loser-winner  | BR-BR |    -0.26 | -0.48, -0.04 | 0.024 
loser-winner  | BR-AR |     0.15 | -0.04,  0.35 | 0.145 
loser-winner  | AR-BE |    -0.06 | -0.09, -0.03 | < .001
loser-winner  | AR-AE |     0.52 |  0.41,  0.63 | < .001
loser-winner  | AR-BR |    -0.36 | -0.47, -0.25 | < .001
loser-winner  | AR-AR |     0.06 |  0.02,  0.09 | < .001
winner-winner | BE-AE |     0.58 |  0.47,  0.68 | < .001
winner-winner | BE-BR |    -0.30 | -0.41, -0.19 | < .001
winner-winner | BE-AR |     0.11 |  0.10,  0.13 | < .001
winner-winner | AE-BR |    -0.87 | -1.02, -0.73 | < .001
winner-winner | AE-AR |    -0.46 | -0.57, -0.36 | < .001
winner-winner | BR-AR |     0.41 |  0.30,  0.52 | < .001
Code
test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

CHECK: Model Assumptions

Code
pp(model);
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
Code
performance::check_model(get(model))

Code
ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)
Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.
Code
ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)
Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.
Code
ls()
 [1] "cat0"                     "check_pkgs"               "control"                  "create_factor_mappings"   "df0"                      "df2"                      "df3"                      "dir8"                    
 [9] "extra"                    "fbase"                    "file"                     "fit01aPh"                 "fit02aPh"                 "fit03aPh"                 "fit04aPh"                 "fpath"                   
[17] "generate_permutations"    "get_stars"                "gg44"                     "gg88"                     "ggpairs_lower_fun"        "hasConverged"             "ifd0"                     "ifn0"                    
[25] "jiko_count_na"            "line0"                    "line1"                    "line2"                    "lineR"                    "model"                    "name"                     "obj_has_name"            
[33] "ofd0"                     "ofd2"                     "pp"                       "pred0"                    "pred0tab"                 "rect0"                    "rect2"                    "rect3"                   
[41] "REML"                     "replace_factor_levels"    "replace_outliers_with_NA" "rescale01"                "rxy0"                     "rxy2"                     "rxy3"                     "sep0"                    
[49] "sep1"                     "sep2"                     "set.w"                    "stuff0"                   "tbl0"                     "terms"                    "test0"                    "test0tab"                
[57] "test2"                    "test2tab"                 "time0"                    "type"                     "var0"                     "write_model_info"        

Auxiliary Checkups

Variables influence / importance in fit04aPh

Code
## TODO CHANGE TITLE
formula(fit04aPh)
Agency ~ (Time | Name) + Time * Phase * Outcome
Code
eff_size0 <- effectsize::eta_squared(fit04aPh)
eff_size2 <- eff_size0 %>%
  dplyr::mutate(Interpret = effectsize::interpret_eta_squared(Eta2_partial))

performance::print_html(eff_size2)
Effect Size for ANOVA
Parameter Eta2 (partial) 95% CI Interpret
Time 2.41e-04 [0.00, 1.00] very small
Phase 6.86e-04 [0.00, 1.00] very small
Outcome 5.68e-03 [0.00, 1.00] very small
Time:Phase 2.04e-03 [0.00, 1.00] very small
Time:Outcome 5.84e-05 [0.00, 1.00] very small
Phase:Outcome 1.01e-03 [0.00, 1.00] very small
Time:Phase:Outcome 1.24e-04 [0.00, 1.00] very small

List models (for Model Hierarchy / Hierarchies)

Code
temp <- lapply(ls(pattern="^fit\\d+"), pp)
fit01aPh: [df0] Agency ~ (1 | Name) + 1
fit02aPh: [df0] Agency ~ (Time | Name) + Time
fit03aPh: [df0] Agency ~ (Time | Name) + Time * Phase
fit04aPh: [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

Compare models using performance package with score

Code
perf0 <- performance::compare_performance(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf0 %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPh lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.49%
fit03aPh lmerModLmerTest 0.11 5.28e-03 0.10 0.26 0.26 4.39e-118 4.39e-118 1.23e-100 51.83%
fit02aPh lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 2.04e-289 2.04e-289 6.95e-259 41.69%
fit01aPh lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.00%
NA
Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- perf0 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  plot = gg88,
  width=8,
  height=6)

gg88

Compare less models

Code
perf2 <- performance::compare_performance(
  ## fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  ## fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf2 %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPh lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 75.00%
fit03aPh lmerModLmerTest 0.11 5.28e-03 0.10 0.26 0.26 4.39e-118 4.39e-118 1.23e-100 25.00%
NA
Code
suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- perf2 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf2.compare-03-score.png"),
  plot = gg44,
  width=8,
  height=6)

gg44  

Comment

Add interpretation of the difference between best fit to data and best predictive power and overfitting (aspecially in the context of the last model).

Code
gg44

Performance Table Sorted by R2_conditional

Code
perf0 %>% dplyr::arrange(desc(R2_conditional)) %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit03aPh lmerModLmerTest 0.11 5.28e-03 0.10 0.26 0.26 4.39e-118 4.39e-118 1.23e-100 51.83%
fit02aPh lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 2.04e-289 2.04e-289 6.95e-259 41.69%
fit04aPh lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.49%
fit01aPh lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.00%
NA

Performance Table Sorted by R2_marginal

Code
perf0 %>% dplyr::arrange(desc(R2_marginal)) %>% performance::print_html()
Comparison of Model Performance Indices
Name Model R2 (cond.) R2 (marg.) ICC RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
fit04aPh lmerModLmerTest 0.10 0.02 0.08 0.26 0.26 1.00 1.00 1.00 84.49%
fit03aPh lmerModLmerTest 0.11 5.28e-03 0.10 0.26 0.26 4.39e-118 4.39e-118 1.23e-100 51.83%
fit02aPh lmerModLmerTest 0.10 4.62e-04 0.10 0.26 0.26 2.04e-289 2.04e-289 6.95e-259 41.69%
fit01aPh lmerModLmerTest 0.08 0.00 0.08 0.26 0.26 0.00e+00 0.00e+00 0.00e+00 1.00%
NA

Interpret R2

Code
model <- "fit01aPh" # [df0] Agency ~ (1 | Name) + 1
model <- "fit02aPh" # [df0] Agency ~ (Time | Name) + Time
model <- "fit03aPh" # [df0] Agency ~ (Time | Name) + Time * Phase
model <- "fit04aPh" # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_conditional, rules="cohen1988"))
weak 
Code
cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_marginal, rules="cohen1988"))
weak 

Plot models: Part 1: Basic Models

Code
suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_models(
  ## CAUTION the null model can not be used here 
  ## Thus to keep the numbers consistent I have 
  ## used model 02 as an input twice
  ## fit01aPh, # [df0] Agency ~ (Time | Name) + Time
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
  spacing=1, dot.size=1
)

ggsave(
  file = file.path(dir8, "summary-all-models-part-01.png"),
  plot = gg88,
  width=5,
  height=4)

gg88

Tab Models

The above `plot_models` seems to be a lot more illustrative

Code
# , eval = TRUE, results='hide'
library(sjPlot)
Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
Code
library(sjmisc)

Attaching package: 'sjmisc'
The following object is masked from 'package:skimr':

    to_long
The following object is masked from 'package:purrr':

    is_empty
The following object is masked from 'package:tidyr':

    replace_na
The following object is masked from 'package:tibble':

    add_case
Code
library(sjlabelled)

Attaching package: 'sjlabelled'
The following object is masked from 'package:forcats':

    as_factor
The following object is masked from 'package:dplyr':

    as_label
The following object is masked from 'package:ggplot2':

    as_label
Code
tab0 <- sjPlot::tab_model(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
    show.reflvl = TRUE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-01-SEE-PNG-images-for-a-better-overview.html"))
tab0
  Agency Agency Agency Agency
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.50 *** 0.49 – 0.50 <0.001 0.50 *** 0.49 – 0.50 <0.001 0.53 *** 0.52 – 0.53 <0.001 0.52 *** 0.51 – 0.52 <0.001
BE Reference Reference Reference Reference
AE 0.02 ** 0.01 – 0.02 0.002 -0.02 * -0.03 – -0.00 0.046
BR -0.13 *** -0.16 – -0.11 <0.001 -0.15 *** -0.20 – -0.09 <0.001
PhaseAE:Outcomewinner 0.05 *** 0.03 – 0.07 <0.001
AR -0.00 -0.01 – 0.01 0.899 -0.10 *** -0.12 – -0.09 <0.001
PhaseAR:Outcomewinner 0.14 *** 0.12 – 0.16 <0.001
loser Reference Reference Reference Reference
winner 0.02 *** 0.01 – 0.04 <0.001
PhaseBR:Outcomewinner 0.02 -0.04 – 0.08 0.472
Time -0.01 *** -0.02 – -0.01 <0.001 0.04 *** 0.04 – 0.05 <0.001 0.06 *** 0.05 – 0.07 <0.001
Time:Outcomewinner -0.02 ** -0.03 – -0.01 0.005
Time:PhaseAE -0.64 *** -0.73 – -0.55 <0.001 -0.88 *** -1.05 – -0.71 <0.001
Time:PhaseAE:Outcomewinner 0.31 ** 0.11 – 0.50 0.002
Time:PhaseAR -0.11 *** -0.12 – -0.10 <0.001 -0.08 *** -0.11 – -0.05 <0.001
Time:PhaseAR:Outcomewinner -0.04 * -0.07 – -0.00 0.033
Time:PhaseBR 0.22 *** 0.13 – 0.32 <0.001 0.02 -0.18 – 0.21 0.860
Time:PhaseBR:Outcomewinner 0.28 * 0.06 – 0.50 0.013
Random Effects
σ2 0.07 0.07 0.07 0.07
τ00 0.01 Name 0.01 Name 0.01 Name 0.00 Name
τ11   0.00 Name.Time 0.00 Name.Time 0.00 Name.Time
ρ01   0.18 Name 0.18 Name -0.06 Name
ICC 0.08 0.10 0.10 0.08
N 870 Name 870 Name 870 Name 870 Name
Observations 169997 169997 169997 169997
Marginal R2 / Conditional R2 0.000 / 0.084 0.000 / 0.102 0.005 / 0.106 0.020 / 0.101
* p<0.05   ** p<0.01   *** p<0.001
Code
## knitr::html
## papaja::
## rmarkdown::render
## rmarkdown::render(tab0$knitr)
## tab0$knitr

Tabulate models (publication ready)

Code
tab0 <- sjPlot::tab_model(
  fit01aPh, # [df0] Agency ~ (1 | Name) + 1
  fit02aPh, # [df0] Agency ~ (Time | Name) + Time
  fit03aPh, # [df0] Agency ~ (Time | Name) + Time * Phase
  fit04aPh, # [df0] Agency ~ (Time | Name) + Time * Phase * Outcome
    show.reflvl = FALSE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-02-publication.html"))
tab0
  Agency Agency Agency Agency
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.50 *** 0.49 – 0.50 <0.001 0.50 *** 0.49 – 0.50 <0.001 0.53 *** 0.52 – 0.53 <0.001 0.52 *** 0.51 – 0.52 <0.001
Time -0.01 *** -0.02 – -0.01 <0.001 0.04 *** 0.04 – 0.05 <0.001 0.06 *** 0.05 – 0.07 <0.001
Phase [AE] 0.02 ** 0.01 – 0.02 0.002 -0.02 * -0.03 – -0.00 0.046
Phase [BR] -0.13 *** -0.16 – -0.11 <0.001 -0.15 *** -0.20 – -0.09 <0.001
Phase [AR] -0.00 -0.01 – 0.01 0.899 -0.10 *** -0.12 – -0.09 <0.001
Time × Phase [AE] -0.64 *** -0.73 – -0.55 <0.001 -0.88 *** -1.05 – -0.71 <0.001
Time × Phase [BR] 0.22 *** 0.13 – 0.32 <0.001 0.02 -0.18 – 0.21 0.860
Time × Phase [AR] -0.11 *** -0.12 – -0.10 <0.001 -0.08 *** -0.11 – -0.05 <0.001
Outcome [winner] 0.02 *** 0.01 – 0.04 <0.001
Time × Outcome [winner] -0.02 ** -0.03 – -0.01 0.005
Phase [AE] × Outcome
[winner]
0.05 *** 0.03 – 0.07 <0.001
Phase [BR] × Outcome
[winner]
0.02 -0.04 – 0.08 0.472
Phase [AR] × Outcome
[winner]
0.14 *** 0.12 – 0.16 <0.001
(Time × Phase [AE]) ×
Outcome [winner]
0.31 ** 0.11 – 0.50 0.002
(Time × Phase [BR]) ×
Outcome [winner]
0.28 * 0.06 – 0.50 0.013
(Time × Phase [AR]) ×
Outcome [winner]
-0.04 * -0.07 – -0.00 0.033
Random Effects
σ2 0.07 0.07 0.07 0.07
τ00 0.01 Name 0.01 Name 0.01 Name 0.00 Name
τ11   0.00 Name.Time 0.00 Name.Time 0.00 Name.Time
ρ01   0.18 Name 0.18 Name -0.06 Name
ICC 0.08 0.10 0.10 0.08
N 870 Name 870 Name 870 Name 870 Name
Observations 169997 169997 169997 169997
Marginal R2 / Conditional R2 0.000 / 0.084 0.000 / 0.102 0.005 / 0.106 0.020 / 0.101
* p<0.05   ** p<0.01   *** p<0.001
Code
# save.image(file=file.path(dir8, "session.RData"))
# load(file.path(dir8, "session.RData")

Report

Code
result <- report::report(fit04aPh)
Code
print(result)
We fitted a linear mixed model (estimated using REML and Nelder-Mead optimizer) to predict Agency with Time, Phase and Outcome (formula: Agency ~ Time * Phase * Outcome). The model included Time as random effects
(formula: ~Time | Name). The model's total explanatory power is weak (conditional R2 = 0.10) and the part related to the fixed effects alone (marginal R2) is of 0.02. The model's intercept, corresponding to Time = 0,
Phase = BE and Outcome = loser, is at 0.52 (95% CI [0.51, 0.52], t(169977) = 115.70, p < .001). Within this model:

  - The effect of Time is statistically significant and positive (beta = 0.06, 95% CI [0.05, 0.07], t(169977) = 10.19, p < .001; Std. beta = 0.12, 95% CI [0.10, 0.14])
  - The effect of Phase [AE] is statistically significant and negative (beta = -0.02, 95% CI [-0.03, -3.46e-04], t(169977) = -2.00, p = 0.046; Std. beta = 0.16, 95% CI [0.06, 0.26])
  - The effect of Phase [BR] is statistically significant and negative (beta = -0.15, 95% CI [-0.20, -0.09], t(169977) = -5.40, p < .001; Std. beta = -0.54, 95% CI [-0.78, -0.30])
  - The effect of Phase [AR] is statistically significant and negative (beta = -0.10, 95% CI [-0.12, -0.09], t(169977) = -10.65, p < .001; Std. beta = -0.37, 95% CI [-0.44, -0.29])
  - The effect of Outcome [winner] is statistically significant and positive (beta = 0.02, 95% CI [0.01, 0.04], t(169977) = 4.08, p < .001; Std. beta = 0.09, 95% CI [0.05, 0.13])
  - The effect of Time × Phase [AE] is statistically significant and negative (beta = -0.88, 95% CI [-1.05, -0.71], t(169977) = -10.21, p < .001; Std. beta = -1.86, 95% CI [-2.22, -1.51])
  - The effect of Time × Phase [BR] is statistically non-significant and positive (beta = 0.02, 95% CI [-0.18, 0.21], t(169977) = 0.18, p = 0.860; Std. beta = 0.04, 95% CI [-0.37, 0.45])
  - The effect of Time × Phase [AR] is statistically significant and negative (beta = -0.08, 95% CI [-0.11, -0.05], t(169977) = -5.24, p < .001; Std. beta = -0.16, 95% CI [-0.22, -0.10])
  - The effect of Time × Outcome [winner] is statistically significant and negative (beta = -0.02, 95% CI [-0.03, -6.25e-03], t(169977) = -2.83, p = 0.005; Std. beta = -0.04, 95% CI [-0.07, -0.01])
  - The effect of Phase [AE] × Outcome [winner] is statistically significant and positive (beta = 0.05, 95% CI [0.03, 0.07], t(169977) = 4.83, p < .001; Std. beta = 0.11, 95% CI [-8.07e-03, 0.23])
  - The effect of Phase [BR] × Outcome [winner] is statistically non-significant and positive (beta = 0.02, 95% CI [-0.04, 0.08], t(169977) = 0.72, p = 0.472; Std. beta = 0.01, 95% CI [-0.27, 0.29])
  - The effect of Phase [AR] × Outcome [winner] is statistically significant and positive (beta = 0.14, 95% CI [0.12, 0.16], t(169977) = 12.52, p < .001; Std. beta = 0.53, 95% CI [0.44, 0.62])
  - The effect of (Time × Phase [AE]) × Outcome [winner] is statistically significant and positive (beta = 0.31, 95% CI [0.11, 0.50], t(169977) = 3.03, p = 0.002; Std. beta = 0.65, 95% CI [0.23, 1.07])
  - The effect of (Time × Phase [BR]) × Outcome [winner] is statistically significant and positive (beta = 0.28, 95% CI [0.06, 0.50], t(169977) = 2.49, p = 0.013; Std. beta = 0.59, 95% CI [0.13, 1.06])
  - The effect of (Time × Phase [AR]) × Outcome [winner] is statistically significant and negative (beta = -0.04, 95% CI [-0.07, -2.93e-03], t(169977) = -2.14, p = 0.033; Std. beta = -0.08, 95% CI [-0.14, -6.19e-03])

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.